@LorenaABarba

Version 0.2 -- March 2014

Flow over a cylinder with source-panels

In previous lessons, we used potential-flow singularities of defined strength to represent the shape of simple geometries, such as a Rankine oval or a circular cylinder, immersed in a free stream. We were rather lucky that when superposing a few fundamental potential-flow solutions, the stream-line pattern that resulted had a closed dividing stream line that we could interpret as a solid body.

But what if we want to represent the stream lines around an arbitrary geometry? Would you be able to define the combination of fundamental solutions to get the expected result? How could you do that? Trial and error? It would take enormous luck and a lot of work to get a geometry we want.

In this lesson, the objective is to calculate the source-strength distribution that can produce potential flow around a given geometry: a circular cylinder. We know that we can get the flow around a cylider by superposing a doublet in a free stream, but here we want to develop a more general approach that can later be extended to different shapes.

The method we will use consists of representing the geometry of the body by a series of short linear segments, called panels, that correspond to source sheets like we studied in the previous lesson.

What we are aiming for is a method that starts with a defined body geometry, then works out the strength of sources needed in each panel representing that geometry to get a dividing streamline right on the body boundary. We will have to solve for the source strengths by specifying that the body be a streamline, i.e., the velocity be tangent there.

Let's start by loading the Python libraries that we will need in this notebook.

We will add a uniform horizontal flow of magnitude u_inf, so let's make that equal to 1:

Definition of the geometry

The geometry considered here will be a circular cylinder of unit radius. We can define this geometry very easily by a set of points going around the angular range between $0$ and $2\pi$.

Discretization into panels

A panel, which represents a source sheet, is defined by its two end-points (xa, ya) and (xb, yb) and its strength sigma. We'll also need its center point (xc, yc) and its length for further calculations. The orientation of the panel is defined by the angle between the $x$-axis and its normal in the counter-clockwise sense.

What information do we need to compute on each panel? First of all, we will need the strength of the source sheet that will lead to the correct streamlines. In addition, we'll also want the tangential velocity (the normal velocity on the body is zero for an inviscid flow) and the pressure coefficient.

In this lesson, you'll really appreciate having learned about classes. It will make the code so much easier to manage. We create a class named Panel containing all the geometry data related to one panel. With a start- and end-point, the class internally calculates the center-point, length and normal vector. It also initializes to zero the source strength, tangential velocity and pressure coefficient. (These will be updated later.)

To store all the discretization, we create a NumPy array of size N_panels where each item in the array is an object of type Panel.

Flow-tangency boundary condition

In Lesson 1, you worked out the velocity potential of a single source as part of the final Challenge Task. It followed simply from integrating the radial velocity, $u_r=\frac{\sigma}{2\pi r}$, that

$$\phi=\frac{\sigma}{2\pi}\ln r$$

(The integration also gives a function of $\theta$ that is seen to be a constant because $u_{\theta}=0$; we take this constant to be zero.)

We will use the velocity potential in this lesson to easily express that the velocity be tangent at the panel, i.e., that $u_n=0$ with:

$$u_n(x,y)=\frac{\partial \phi}{\partial n}(x,y)$$

at a given point on the panel. We choose the point to enforce that velocity-tangency condition as the center of the panel (and we call it the control point).

The velocity potential in Cartesian coordinates of a source sheet on a panel is

$$\phi\left(x,y\right) = \frac{\sigma}{2\pi} \int_\text{panel} \ln \sqrt{ \left(x-x(s)\right)^2 + \left(y-y(s)\right)^2 } {\rm d}s$$

where $s$ is the running coordinate along the panel and $\left(x(s),y(s)\right)$ are the Cartesian coordinates of $s$.

Superposition of the potential of each panel gives the total potential at any point $\left(x,y\right)$, so we make a sum of all the panel contributions as follows (moving the $\frac{1}{2}$ exponent in the logarithmic term as a factor outside the integral):

$$\phi\left(x,y\right) = \sum_{j=1}^{N_p} \frac{\sigma_j}{4\pi} \int \ln \left( \left(x-x_j(s_j)\right)^2 + \left(y-y_j(s_j)\right)^2 \right) {\rm d}s_j$$

By finally superposing the free stream, the flow around an immersed circular cylinder will be represented by the following velocity potential:

$$\phi\left(x,y\right) = U_\infty x + \sum_{j=1}^{N_p} \frac{\sigma_j}{4\pi} \int \ln \left( \left(x-x_j(s_j)\right)^2 + \left(y-y_j(s_j)\right)^2 \right) {\rm d}s_j$$

Enforcing the flow-tangency condition on each control point approximately makes the body geometry correspond to a dividing streamline (and the approximation improves if we represented the body with more and more panels). So, for each panel $i$, we make $u_n=0$ at $(x_{c_i},y_{c_i})$:

$$u_{n_i} = \frac{\partial}{\partial n_i}\left\lbrace \phi\left(x_{c_i},y_{c_i}\right) \right\rbrace = 0$$

which leads to

$$0 = U_\infty \cos\beta_i + \sum_{j=1}^{N_p} \frac{\sigma_j}{2\pi} \int \frac{ \left(x_{c_i}-x_j(s_j)\right) \frac{\partial x_{c_i}}{\partial n_i} + \left(y_{c_i}-y_j(s_j)\right) \frac{\partial y_{c_i}}{\partial n_i} } {\left(x_{c_i}-x_j(s)\right)^2 + \left(y_{c_i}-y_j(s)\right)^2} {\rm d}s_j$$

where $\beta_i$ is the angle that the panel's normal makes with the $x$-axis, so

$$\frac{\partial x_{c_i}}{\partial n_i} = \cos\beta_i \quad \text{and}\quad\frac{\partial y_{c_i}}{\partial n_i} = \sin\beta_i$$

and

$$x_j(s_j) = x_{a_j} - \sin\left(\beta_j\right) s_j$$$$y_j(s_j) = y_{a_j} + \cos\left(\beta_j\right) s_j$$

But, there is still a problem to handle when $i=j$. From the previous notebook, we have seen that the stength of the source sheet should be a specific value so that the streamlines do not penetrate the panel. This helps us determine that the contribution of the $i$-th panel to itself is $\frac{\sigma_i}{2}$.

Finally, the boundary condition at the center point of the $i$-th panel gives

$$0 = U_\infty \cos\beta_i + \frac{\sigma_i}{2} + \sum_{j=1,j\neq i}^{N_p} \frac{\sigma_j}{2\pi} \int \frac{ \left(x_{c_i}-x_j(s_j)\right) \cos\beta_i + \left(y_{c_i}-y_j(s_j)\right) \sin\beta_i } {\left(x_{c_i}-x_j(s)\right)^2 + \left(y_{c_i}-y_j(s)\right)^2} {\rm d}s_j$$

From the equation above, we understand that we will have to compute integrals using the SciPy function integrate.quad(). We define a function integral_normal() that will do the job.

Solving the system of equations

We just developed an equation to enforce a flow-tangency condition on the $i$-th panel. There are N_panels panels $i$ and N_panels unknown strengths $\sigma_i$. Therefore, the problem represents solving a linear system of equations of the form

$$[A][\sigma] = [b]$$

where

$$A_{ij} = \begin{cases} \begin{matrix} \frac{1}{2} & \mbox{, if } i=j \cr \frac{1}{2\pi} \int \frac{ \left(x_{c_i}-x_j(s_j)\right) \cos\beta_i + \left(y_{c_i}-y_j(s_j)\right) \sin\beta_i } {\left(x_{c_i}-x_j(s)\right)^2 + \left(y_{c_i}-y_j(s)\right)^2} ds_j & \mbox{, if } i\neq j \end{matrix} \end{cases}$$

and

$$b_i = - U_\infty \cos\beta_i$$

for $1\leq i,j \leq N_p$. Let's fill a matrix $A$ and a right-hand-side vector $b$ with the necessary values:

Hey! We just used a new Python built-in function: enumerate(). It allows us to have access to each element panel in the array panels while keeping a count i (that starts from 0) to locate the element of A to fill.

Now, we can easily solve the linear system of equations using the function linalg.solve() from NumPy, and assign each source-panel its appropriate strength:

Pressure coefficient on the surface

At this point, we have the source strength distribution required to compute the streamlines around our geometry. A very useful measurement of the results is the pressure coefficient along the surface of the geometry.

From Bernoulli's equation, the pressure coefficient on the $i$-th panel is

$$C_{p_i} = 1-\left(\frac{u_{t_i}}{U_\infty}\right)^2$$

where $u_{t_i}$ is the tangential component of the velocity at the center point of the $i$-th panel,

$$u_{t_i} = \frac{\partial}{\partial t_i}\left\lbrace \phi\left(x_{c_i},y_{c_i}\right) \right\rbrace$$

which we can obtain as:

$$u_{t_i} = -U_\infty \sin\beta_i + \sum_{j=1}^{N_p} \frac{\sigma_j}{2\pi} \int \frac{ \left(x_{c_i}-x_j(s_j)\right) \frac{\partial x_{c_i}}{\partial t_i} + \left(y_{c_i}-y_j(s_j)\right) \frac{\partial y_{c_i}}{\partial t_i} } {\left(x_{c_i}-x_j(s)\right)^2 + \left(y_{c_i}-y_j(s)\right)^2} {\rm d}s_j$$

with

$$\frac{\partial x_{c_i}}{\partial t_i} = -\sin\beta_i \quad\text{and} \quad \frac{\partial y_{c_i}}{\partial t_i} = \cos\beta_i$$

Note that the contribution to the tangential velocity at a source panel from its own velocity potential is zero, which makes sense because streamlines go outwards from a source.

We define a function integral_tangential() that will compute the integrals above using the SciPy function integrate.quad() once again:

Once we have computed the tangential velocity on each panel, we can calculate the pressure coefficient.

Allright! It is time to plot the surface pressure coefficient.

Just before this, we should remember that in the lesson on the doublet, we found that the exact pressure coeficient on the surface on a cylinder was

$$Cp = 1 - 4\sin^2 \theta$$

i.e.

$$Cp = 1 - 4\left(\frac{y}{R}\right)^2$$

We can use this to compare with the results obtained with our source-panel code.

Challenge task

Now that we have computed the pressure coefficient on the surface of the cylinder, it will be interesting to visualize what the streamlines look like.

To do that, we use the function streamplot() from Matplotlib, requiring the Cartesian velocity components (u,v) on a mesh grid (X,Y). Therefore, the first step is to derive the equations for the velocity components.

The potential at point $\left(x,y\right)$ of the $N_p$ source sheets in a uniform horizontal flow $U_\infty$ is

$$\phi\left(x,y\right) = U_\infty x + \sum_{j=1}^{N_p} \frac{\sigma_j}{4\pi} \int \ln \left( \left(x-x_j(s_j)\right)^2 + \left(y-y_j(s_j)\right)^2 \right) {\rm d}s_j$$

And the velocity field at point $\left(x,y\right)$ is

$$u\left(x,y\right) = \frac{\partial}{\partial x}\left\lbrace \phi\left(x,y\right) \right\rbrace$$$$v\left(x,y\right) = \frac{\partial}{\partial y}\left\lbrace \phi\left(x,y\right) \right\rbrace$$

Your task will be to:

• derive the Cartesian velocity components
• create a mesh grid
• compute the velocity field on the mesh grid
• plot the results
• change the number of panels to improve the visualization

---